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RAPID GENERATION OF OPTIMAL ASTEROID POWERED 
DESCENT TRAJECTORIES VIA CONVEX OPTIMIZATION 

Robin Pinson* and Ping Lu f 

This paper investigates a convex optimization based method that can rapidly generate the fuel 
optimal asteroid powered descent trajectory. The ultimate goal is to autonomously design the 
optimal powered descent trajectory on-board the spacecraft immediately prior to the descent 
burn. Compared to a planetary powered landing problem, the major difficulty is the complex 
gravity field near the surface of an asteroid that cannot be approximated by a constant gravity 
field. This paper uses relaxation techniques and a successive solution process that seeks the 
solution to the original nonlinear, nonconvex problem through the solutions to a sequence of 
convex optimal control problems. 

INTRODUCTION 

Mission proposals that land spacecraft on asteroids are becoming increasingly popular. However, in order 
to have a successful mission the spacecraft must reliably and softly land at the intended landing site with pin- 
point precision. The problem under investigation is how to design a fuel optimal powered descent trajectory 
that can be quickly computed on-board the spacecraft, without interaction from ground control. An optimal 
trajectory designed immediately prior to the powered descent bum has many advantages. These advantages 
include the ability to use the actual vehicle starting state as the initial condition in the trajectory design and the 
ease of updating the target landing site. For long trajectories, the trajectory can be updated periodically by a 
redesign of the optimal trajectory based on current vehicle conditions to improve guidance performance. One 
of the key drivers for being completely autonomous is the infrequent and delayed communication between 
ground control and the vehicle. Challenges that arise from designing an asteroid powered descent trajectory 
include complicated nonlinear gravity fields, small rotating bodies and low thmst vehicles. There are many 
factors that will not be understood until the spacecraft reaches the asteroid, which is why many missions 
spend long periods of time near the asteroid characterizing it and choosing the landing site. It is imperative 
to make the trajectory design algorithm as flexible as possible to account for this new information. 

There have been two successful landings on asteroids and one on a comet over the last two decades. NEAR 
was the first successful landing. NASA decided to attempt a controlled descent, after debating options for the 
spacecraft following completion of the primary mission objectives. The NEAR spacecraft was not designed to 
land on Eros, which led to designing landing operations after launch. Four braking maneuvers were planned 
over a 50 minute period. Each maneuver targeted a specified decrease in vehicle velocity. It was expected that 
the vehicle would impact the surface with a 1.3 m/s velocity. The spacecraft impacted earlier than expected 
in the range of 1.5 -1.8 m/s. 1 Due to malfunctioning equipment, Hayabusa was painstakingly guided to a 
point near the landing site from ground control, even with a 20 minute communication delay. Close to the 
surface, landing target markers were dropped to the surface and Hayabusa rendezvoused with the markers to 
successfully reach the surface. 2 Rosetta’s lander Philae followed a ballistic trajectory down to the comet’s 
surface. When the trajectory was originally analyzed there were two options, a completely passive trajectory 
and an option containing one pre-planned maneuver. 3,4 

There are a handful of proposed missions to asteroids whose landing algorithms and missions are still 
being finalized. Most rely on ballistic trajectories or ground designed trajectories and waypoints. The lander 
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MASCOT relies on the primary spacecraft to drop down to 100 m altitude and then MASCOT free falls to the 
surface for the Marco Polo and Hayabusa-2 missions. 5,6 OSIRIS-REx has a checkpoint and a matchpoint at 
which maneuvers will be performed. These are targeted maneuvers which are adjusted on-board to account 
for differences between the originally designed points and the actual spacecraft conditions at those points. 
The last maneuver at the matchpoint is eight minutes prior to touchdown at 45 m altitude and slows the 
descent velocity. 7 

Since the asteroid powered descent trajectory design is a new field of study, there are a wide range of 
proposed methods and ideas that are being developed. All of these ideas have advantages and disadvantages. 
The trajectory design and related guidance algorithms are both studied, as some methods forgo trajectory 
design and focus entirely on closed loop guidance efforts. One proposed approach to the rapid generation 
of the descent trajectory is a semi-analytical approach which assumes the acceleration profile is a cubic 
polynomial. Given the acceleration profile shape with known initial and final points, the problem is solved 
analytically except for three parameters: time of flight, initial thrust and initial thrust angle. Combining 
these with path constraints yields a nonlinear program which is then solved. This method does allow for 
rapidly designing a trajectory at the cost of a non-optimal solution due to the acceleration assumption. 8 A 
second approach for designing the descent trajectory focuses on extracting information from the optimal 
control theory formulation of the problem. An indirect method that propagates and solves the costates, 
boundary conditions and differential equations is combined with a direct method that optimizes the original 
cost function. This combination decreases the computation time as compared to directly optimizing the 
nonlinear cost function and original problem formulation. 9 

Variations on sliding mode control theory as a guidance algorithm for the descent trajectory is found in 
several places. Sliding mode control theory is a robust control law when faced with bounded disturbances. 
One method uses multiple sliding surfaces to drive the spacecraft to the landing site in a finite amount of 
time. 10 Other proposals use a predesigned trajectory and apply variations on sliding mode control theory to 
counteract disturbances while tracking the trajectory. 11,12 One proposed method of predesigning the trajec- 
tory is to apply a homotopic and successive solution technique while solving the state and costate differential 
equations via a shootout method. 12 

There are two previous studies that form the background to the current investigation. The first study 
looked in-depth at applying convex optimization to a powered descent trajectory on Mars with promising 
results. 13 This showed that the powered descent equations of motion can be relaxed and formed into a convex 
optimization problem and that the optimal solution of the relaxed problem is indeed a feasible solution to 
the original problem. This analysis used a constant gravity field, which is a poor approximation for an 
asteroid. In addition, these authors have looked into including thrust pointing constraints, while achieving 
lossless convexification. 14 The second study applied a successive solution process to formulate a second order 
cone problem that designs rendezvous and proximity operations trajectories. 15 These trajectories included a 
Newtonian gravity model. The equivalence of the solutions between the relaxed and the original problem was 
theoretically established. 

Convex optimization had great success in the above studies and the focus of this analysis is the applicability 
to the asteroid powered descent problem. A key advantage of using convex optimization and its subclasses 
is that as long as there is at least one feasible solution, the global minimum is guaranteed in a finite number 
of steps. A wide range of functions and formulations are included in the convex variable class and convex 
optimization including second order cone programming and linear programming. 

The fuel optimal control problem in this work is to determine the optimal (finite) thrust vector (both magni- 
tude and direction) to land the spacecraft at a specified location, in the presence of a highly nonlinear gravity 
field and subject to various mission and operational constraints. The proposed solution is to use convex opti- 
mization, a gravity model with higher fidelity than Newtonian, and an iterative solution process to design the 
fuel optimal trajectory. The solution to the convex optimization problem is the thrust profile, magnitude and 
direction, that will yield the minimum fuel trajectory for a soft landing at the target site. 

This paper discusses the optimization problem formulation and spacecraft dynamics, including the gravity 
model. The original formulation is not in a convex optimization format. The steps and methods taken to 


2 



generate the problem into a convex form that will yield the same optimal solution as the original problem 
is the main focus. Relaxation techniques, investigation into problem equivalence, change of variables and a 
successive solution method are all used to generate the final optimization problem that can be solved with a 
convex optimizer. Following this discussion will be results and trends seen from generating trajectories with 
this process. 

ORIGINAL PROBLEM FORMULATION 

For a soft landing on an asteroid, the trajectory is designed as a two point boundary value problem. The 
state of the spacecraft at the start of the powered descent burn is known, along with the landing site location. 
A soft landing occurs when the spacecraft has zero velocity relative to the landing site when it reaches the 
landing site. 



Figure 1. Asteroid centered fixed coordinate system example. Scale is in meters. 


The dynamics of the landing trajectory and the corresponding equations of motion are formulated in an 
asteroid centered fixed (ACF) Cartesian coordinate system. This coordinate system, depicted in Figure 1, 
is fixed at the asteroid’s center of mass. The x axis is aligned with the largest semi-major axis (minimum 
moment of inertia), the z axis is aligned with the smallest semi-major axis (maximum moment of inertia) and 
the y axis is aligned with the intermediate semi-major axis. For asteroids that are not perfect ellipsoids, the 
axes are as closely aligned as possible, while still keeping an orthogonal coordinate system. The majority of 
the work will be in the Cartesian coordinate system; however, some models are originally derived in spherical 
coordinates and then transformed into Cartesian coordinates. The variable r is the radius from the center of 
the asteroid to the spacecraft. The two angles are latitude ((5) measured from the z axis to the radius vector 
and longitude (A) measured from the x axis to the projection of the radius in the x-y plane. 


The ACF is a rotating frame around the asteroid’s spin axis cc with a spin rate of u). In order to apply 
the dynamic equations to all asteroids, including those with a rotation axis not aligned with the z axis or 
tumbling asteroids, the spin axis is allowed to be a time-varying vector. The equations of motion of the 
spacecraft include the rotational effects of the coordinate system, Eq. (1). 


T + m\7U = rn 


( dv _ _ . 

[—+ 2oj x v + oj x r + oj x (c d x r) 


( 1 ) 


The position and velocity vectors of the spacecraft in the ACF frame are represented by r and v. The vehicle 
thrust vector is T and the acceleration due to gravity is the gradient of the gravitational potential (VC/). 


The gravitational potential ( U ) is a function of spacecraft position for all the gravity models: Newtonian, 
spherical potential, polyhedron and interior spherical potential. For this paper, the spherical potential function 
for a triaxial ellipsoid will be used as the gravitational potential. The standard form of the spherical potential 
is given in Eq. (2). 16 

00 ^ i 

u = (~) Pl ’ m [ sin £] {Ci,m cos (mA) + Si t m sin (toA)} (2) 

1=0 m = 0 
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In the model located in Eq. (2), the gravitational potential is based on the distance from the center of mass 
or radius (r), the latitude (5) and the longitude (A) of the spacecraft. The coefficients C) ?m and are 
unique to the asteroid. These coefficients can be determined in a number of ways. The reference radius (r 0 ) 
is associated with the coefficients, usually taken to be the maximum radius of the asteroid. The gravitational 
constant is represented by n and P^ m is the associated Legendre function. While the potential function is 
an infinite sum in order to represent the effects of the mass over the entire body shape, the summation can 
be truncated and still yield reasonable accuracy. The desired accuracy dictates the order and degree of the 
model that is used. A good initial model for an asteroid is to use a 2x2 model, and assume the asteroid is a 
homogenous (constant density) triaxial ellipsoid. 17 For this case, the upper limits on the summation are both 
2. Due to the fact that the origin of the ACF coordinate system is at the center of mass and the symmetry of 
the homogenous ellipsoid, only Co,o, C'2,0 and 62,2 are nonzero. By definition Co,o = 1. The gravitational 
potential for the triaxial ellipsoid reduces to Eq. (3). 


u = ± 


1 + (— Y 

V r / 


C'2,0 ( ^ sin 2 (5 - 1 ) + 6*2,2 (3 cos 2 < 5 ) cos (2A) 


]} 


(3) 


The acceleration components due to gravity or the gradient of the gravitational potential in the ACF Carte- 
sian coordinate system are given in Eq. (4). 
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The terms G i, G 2 and G3 are defined in Eq. (5). 


r 2 

G\ =-| {3G 2 ,o [l.5,sm 2 (6) - 0.5] + 9 C 2 , 2 cos 2 (5) cos (2A)} 
r 2 

G 2 {1.5(72, osin (28) - 3C 2 ,2sin (28) cos (2A)} 

G 3 =f 2 | ( r 2 + ~ 2 ) 6C 2,2 cos 2 (5) Sin (2A) | 


(4) 


(5) 


The gradient in Eq. (4) is highly nonlinear in the position vector, including radius raised to the fifth power 
and cross products of position components. When C'2,0 and 62,2 are zero, Eq. (4) reduces to the standard 
Newtonian field with ^ = — -j^r. 

The problem’s optimization objective is minimization of the fuel usage (or propellant usage), while still 
landing softly. This is achieved by maximizing the mass of the vehicle when it touches down, thus minimizing 
fuel consumption for the given initial mass. To accomplish this, the minimization cost function (Eq. (6)) is the 
negative of the mass at the landing site. The full optimization problem is stated in Eq. (6) through Eq. (18), 
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referred to as Problem PI. The states are f, v, m and the control vector is T . 


min —m(tf) 
s.t. r= v 


( 6 ) 

(7) 


v = 2uj xv — ujxf— ujx (w x r ) + VU (r) 
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( 10 ) 


||r|| cos# — r T 77 < 0 (11) 

777- ^ 777 dry (12) 

f (0) = f 0 (13) 

(' (0) = v 0 (14) 

777 (0) = m we t (15) 

r(tf) = f f (16) 

v(tf)=v f (17) 

tf given (18) 


Initial conditions, position (r 0 ), velocity (f7 0 ) and total spacecraft mass including propellant (m wet ), are 
known, Eq. (13), Eq. (14) and Eq. (15). This is a fixed final time (tf) problem with tf specified. The 
final position (ff) and velocity (ff) are based on the landing site, Eq. (16) and Eq. (17). In the problem 
formulation, once the engine has been turned on, it is assumed to be throttleable between a minimum and 
maximum magnitude (T m ^ n and T ma;r ) with T m ^ n > 0, because the engine cannot be turned off until the 
descent burn is complete. This is accounted for by inequality constraints on the thrust magnitude, Eq. (10). 
The associated exit velocity (v ex ) is the product of the vehicle Isp and reference gravity at Earth sea level. The 
final mass is free, as long as the vehicle does not use more than the allowable propellant which is protected 
for by the constraint in Eq. (12). When the vehicle mass reaches the vehicle dry mass, all propellant 

has been consumed. 

The above optimal control problem is nonlinear and nonconvex. The presence of the thrust vector norm in 
Eq. (9) and Eq. (10) and the mass in the denominator of Eq. (8) cause those equations to be nonlinear. The 
gravitational acceleration in Eq. (8) is highly nonlinear in the position vector. The thrust constraint Eq. (10) 
can be divided into two constraints. 


< T 

— L max 


(19) 


The second inequality in Eq. (19) is convex and an allowable inequality for a convex optimization problem, 
while the first inequality is not. Therefore, to solve this problem by convex optimization, which requires 
equality constraints to be linear and inequality constraints to be convex, the nonlinearities in this problem 
will need to be handled with a variety of techniques. 
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The constraint in Eq. (11) is a glide slope constraint, which represents an approach cone surrounding the 
landing site. The cone is defined by the angle 6 , which is limited to 0 < 6 < 90 deg. When 6 is 90 deg, 
this ensures that the vehicle altitude is always positive and the trajectory design does not try to fly below the 
surface before reaching the landing site. The unit vector normal to the asteroid surface at the landing site is 
n. The angle between the spacecraft position vector and the landing site normal is a , which the constraint 
ensures to be less than or equal to 0. 


CONVEXIFICATION AND RELAXATION OF THE PROBLEM 


The magnitude of the thrust vector or norm of the thrust vector is one of the nonlinearities in the original 
problem. In order to remove nonlinearities, a slack variable T m will be added to become a fourth control to 
supplement T. In the relaxed optimal control problem, this variable is completely separate from the thrust 
vector and only related through a new constraint, Eq. (20). 



<T r 


( 20 ) 


This slack variable T m replaces the magnitude of the thrust vector in Problem PI. The new optimal control 
problem, now called Problem P2, is defined in Eq. (21) through Eq. (33). 


min — mn(tf ) 

(21) 
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(27) 

f (0) = ro 

(28) 

II 

(29) 

m (0) = m wet 

(30) 

r(t f ) = r f 

(31) 

v(tf) = v f 

(32) 

tf given 

(33) 


With the introduction of the slack variable T m Eq. (24) and Eq. (26) are now linear and the constraint in Eq. 
(25) is convex. The remaining nonlinearities are the first and last terms on the right hand side of Eq. (23). 
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The glide slope constraint in Eq. (11) is not included in Problem P2 for brevity of the discussion below. 
Correspondingly, it is understood that the glide slope constraint has also been removed from Problem PI in 
the following discussion. 


As Problem P2 will be solved in this paper, an immediate question arises as to the relationship between the 
solutions to Problems PI and P2. To answer this critical question, examine the following results. 


Proposition 1: Problems PI and P2 are identical when the inequality in Eq. (25) is active, i.e., 
Moreover, the solution to Problem P2 is the same as the solution to Problem PI. 


f 


= Try 


Proof: The first part of Proposition 1 is straightforward. It is immediately clear that when the constraint in 
Eq. (27) is active, Problem P2 reduces to Problem PI. To prove the second part of the Proposition, apply the 
Minimum Principle from optimal control theory 18 and form the Hamiltonian in a vector and matrix format 
that removes the cross products, Eq. (34). Note that the cost function is only dependent on the mass at the 
final time, </>(£/) = The dynamics of the system are autonomous causing the Hamiltonian to be 

constant during the entire flight. 
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(34) 


The states are r, u, m, and the corresponding costates are p r ,Vv, Pm • The control vector ^ consists of the 
thrust vector and slack variable and the constraints on the control vector are described in Eq. (35). 


= 


{(f,r ro ) 


<C P P < p p < p 

^ i m) ± min _ rn 5 J -m _ - 1 rr 


(35) 


From optimal control theory, the differential equations for the costates are listed in Eq. (36) through Eq. (38). 

Pr = W?p v + Wlp v - pj V 2 f/(r) (36) 

Pv = -Pr + 2 wf Pv (37) 


Pm =PyT ( 2 - 
m 2 


(38) 


The problem is a fixed final time problem. Applying the end point constraints and the transversality conditions 
yields one useful equation, Eq. (39). 

Pm (tf) = - 1 (39) 


From the Minimum Principle, the optimal solution is determined by minimizing the Hamiltonian with 
respect to the control vector (^ ) subject to the constraints in Eq. (35). The Hamiltonian is linear in terms of 
the control vector, therefore the minimum will be on the boundary of the admissible control set. This can be 
turned into a constrained optimization problem which minimizes the Hamiltonian by applying the Kasrush- 
Kuhn-Tucker (KKT) conditions. First, the terms in the Hamiltonian expression that are not dependent on 
the control vector can be removed from the optimization problem, as they have no effect on the solution. 
These terms are pjv + (— 2W\v — Wfr — W$r + VC/ (r)). It is interesting to note that as long as the 

gravitational potential remains a function of spacecraft position, it does not influence the minimization of the 
Hamiltonian. Increasing the complexity of the gravity model will not invalidate these arguments. 

The point- wise minimization problem with respect to T, P m is stated in Eq. (40). 


s.t. 
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min ±pj? - ^p m T m 

< P P < P P < P 
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(40) 
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This is a convex optimization problem which has a global optimal solution. The Lagrangian for the con- 
strained optimization problem is listed in Eq. (41). 

— — Pv T PmTm + Ai ( T — Tm) + A2 {T m i n — T m ) + A3 (T m — T rnax ) ( 41 ) 

m v ex V / 

The variables Ai, A2, A3 are the Lagrange multipliers associated with the inequality constraints on T, T m . 
From the KKT conditions, a constraint is active when A^ > 0, while a constraint can only be inactive when 
A i = 0. Applying the KKT conditions yields Eq. (42) and Eq. (43). 


9Sf 

dT 

d jgf 


— — Pv + A 
m 



— Pm ^1 

Vex 


= 0 

— A2 + A3 
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(42) 

(43) 


Now examine the problem as two specific cases when p v is nonzero almost everywhere (a.e.) in [0, tf] and 
when it is zero in a finite time interval. For the first case it will be shown that Ai 7 ^ 0, therefore the solutions 
to Problems PI and P2 are the same. Then by contradiction, show that p v cannot be zero in a finite time 
interval. 

Suppose that p v / 0 in any finite time interval. From Eq. (42), when T / 0, then Ai must be nonzero 
(Ai 7 ^ 0); otherwise, Eq. (42) necessitates p v = 0, since mass is nonzero, contradicting the assumption 
that p v 7 ^ 0. When T = 0,T/||T||isa finite vector. The same argument applies to ascertain that Ai 7 ^ 0. 

Therefore, the corresponding constraint must be active, T = T m . Consequently, the solution to Problem 
P2 is also the solution to Problem PI. 

Next look at the case when p v = 0 for a finite time period. If p v = 0 in a finite time interval, its derivative 
must also be zero simultaneously. From the costate equations in Eq. (36) through Eq. (38), this implies that 
p r = 0 and Pm is a constant. Because the costate system in Eq. (36) through Eq. (38) is a homogeneous 
linear system in p r and p v , it follows that p r = p v = 0 in the entire interval [0, tf] , and p m is constant. From 
the transversality condition in Eq. (39), p m (t) = — 1. The Hamiltonian reduces to Eq. (44). 

H = —T m (44) 

Vex 

The minimum of the Hamiltonian with respect to is the smallest admissible T m , which is T m in > 0, 
therefore T m = Tmin- Thus the Hamiltonian becomes Eq. (45). 

H = T Tmin (45) 

Vex 

If p v = 0 in a finite interval, the constancy condition of the Hamiltonian, Eq. (45), should always hold, 
regardless of the initial conditions, landing site, or the specified final time tf. In particular, if the value of 
tf happens to be equal to tp the optimal flight time in a free-time problem, which is the same as Problem 
P2 except that the final time is free, the Hamiltonian is zero along the optimal trajectory from the Minimum 
Principle. The nonzero constancy condition in Eq. (45) then constitutes a contradiction in this case. There- 
fore, the condition that p v = 0 in a finite interval cannot happen. Since this case is invalid, the case of p v / 0 
a.e. on [0, tf] is the only possibility, which in turn means that T = T m in the solution to Problem P2. This 
solution is also the solution to Problem PI. 


CHANGE OF VARIABLES 

The next step of the process that turns Problems PI and P2 into a convex optimization problem is a change 
of variable for two variable types. The first changes the thrust variables into acceleration due to thrust with 
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the conversion seen in Eq. (46). 


a t = 

atm 


T 

m 

Trn 

m 


(46) 


In addition, a new mass variable and its derivative are introduced to replace the current mass variables. This 
transformation is seen in Eq. (47). 


q = In (rri) 

. 

Vex 


(47) 


This change affects the cost function. The properties of the natural log function allow maximum mass to 
correspond to the maximum of q. 

After the change of variables, the constraint in Eq. (26) becomes the two inequality expressions in Eq. 
(48). 


T-'min^’ ^ — afm 
atm — Tmax^ ^ 


(48) 


Taylor’s series expansion will be applied to e~ q , in order to turn the top inequality into a quadratic equation 
in state and the bottom into a linear. This effectively creates a second order cone constraint and a linear 
constraint, which are both convex. A quadratic cannot be used on the bottom constraint as that would not 
yield a convex constraint. The new constraints (Eq. (49)) are equivalent to the old as long as the expansion 
does not deviate from the original mass. Spot checking shows the difference is on the order of a few grams, 
which is negligible to ensuring the thrust bounds are not exceeded. 


1 - (q ~ Qo) + 0.5 (q - q 0 Y 

atm — Tmax^ ^ ° [1 (O. Qo)\ 


^ atm 


(49) 


The point about which the expansion occurs (q Q ) is calculated ahead of time. It is vehicle mass profile 
assuming the vehicle has been running at minimum thrust for the entire flight, Eq. (50). 


i I J -mm A , 

q Q = In m wet At 


(50) 


In the above equation, At is the time elapsed from the start of the bum. 

After the change in variables and Taylor series expansion, the optimization problem becomes Problem P3 
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found in Eq. (51) through Eq. (64). 


min— q(tf) (51) 

s.t. f=v (52) 

v = a t — 2Co x v — Co x r — Co x (Co x r) + VU (r) (53) 

Q = ~— atm (54) 

Vex 

\\^t || tttm (55) 

T m i n e~ qo 1 - (q ~ q 0 ) + 0.5 (q - q Q ) 2 < a tm (56) 

atm < T max e~ q ° [1 - (q - q Q )] (57) 

q > Qdry (58) 

f (0) = r 0 (59) 

V (0) = v 0 (60) 

q (0) = In ( m wet ) (61) 

r(tf)= r f (62) 

v ( tf ) = Vf (63) 

t f given (64) 


Problem P3 is equivalent to Problem P2 as long as the Taylor series expansion is a good approximation. 
The optimal solution for Problem P3 is the optimal solution of Problem P2 and thus Problem PI. With the 
exception of the gravitational acceleration, the dynamics are linear in state and control, while the inequality 
constraints are linear or second order cones, yielding a convex optimization problem. The gravitational 
acceleration nonlinear term will be handled via a successive solution method. 


SUCCESSIVE SOLUTION METHOD 

The gravitational acceleration is the most challenging aspect to manipulate into a convex form. Recall 
that the 2x2 triaxial ellipsoid gravitational acceleration is highly nonlinear in the position vector including 
radius raised to the fifth power and cross multiplication of position terms. The best way to handle these 
nonlinearities is to introduce a successive solution method. In the successive solution method a series of 
convex optimization problems are solved with each one using data from the previous solution, also known as 
iterations. The equations of motion from Problem P3, found in Eq. (52) through Eq. (54), can be rearranged 
into Eq. (65), where (k) is the current iteration and (k-1) is the previous iteration. 


£(fc) = A ( r ( fc -i)) x (k) + Bu (k) + c (r^)) (65) 

The full state is represented by x = [r T , v T , q] T and the control vector u = [a t T , a trn \ T . The B matrix 
is constant throughout the flight and between iterations. The terms from the gravitational acceleration are 
divided between the A matrix and c vector and evaluated with the previous iteration’s solution. Five different 
arrangements of the terms in A and c were tried. The best one placed the Newtonian terms in the A matrix 
and all the remaining terms in the c vector, see Eq. (66). Thus, c can be thought of as a vector of higher 
order terms. With this arrangement it should be simple to swap out the 2x2 triaxial model for a higher fidelity 
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model. 
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A and c are functions of position, since the nonlinear terms in the gravitational acceleration are position. 
For the first iteration the position vector is taken to be the initial condition for the entire trajectory, equivalent 
to assuming the spacecraft hovered at the initial point for the entire length of trajectory. Iterations continue 
until the current iteration and previous iteration match within a specified tolerance. At this point they are 
close enough to be called equivalent and the gravitational acceleration used to develop the current trajectory 
is no longer an approximation. 

By using the previous trajectory to evaluate the gravitational acceleration, the dynamics are now linear and 
thus the problem is a convex optimization problem. Actually, it is in a second order cone program format, 
which is a subset of convex optimization. In order to keep the problem generic and allow for additional 
constraints or dynamics, a convex optimizer will be used. A second order cone program is also an option for 
solving the problem in its current formulation. 


SCALING AND DISCRETIZATION 

Scaling is important to ensuring good and consistent convergence in all cases. The dynamics of the sys- 
tem are scaled to keep the values near one and to keep any variable from having extra emphasis due to its 
magnitude. Originally, the problem was run without scaling and the optimizer had difficulties finding a fea- 
sible solution. All distance parameters are scaled by Rq, velocity by v sc , acceleration by go, time by t sc and 
angular rate by c o sc . These scale factors are seen in Eq. (67). 

#0 = 7, 90 = ^ 2 , v sc = \/-Ro5o , tsc*.yp^, u sc =^j^ (67) 

The distance scale factor is set to the smallest semi-major axis value (7). This choice is consequential, 
because it keeps the scaled radius to no smaller than unity, which in turn will prevent the terms such as 1 /r 3 
and 1/r 5 in the gravitational acceleration from becoming very large values. When the largest semi-major 
axis was used, the iterations failed to converge on a single answer. Instead the iterations bounced between 
two trajectories. The mass variable (q) is not scaled as the natural log function acts as a scaling function and 
keeps q near one. 

The linear equations of motion are in continuous time; however, optimizers require the problem to be in 
discrete time. A simple discretization and propagation with the trapezoidal rule was applied. The problem 
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is divided into n number of time points with a fixed dt time step between each point. The trapezoidal rule is 
used to propagate the trajectory between the points as demonstrated in Eq. (68). 

Xj — Xj—\ T" 0.5 dt (AjXj H - Aj—^Xj—i H - BjUj A B j A Cj H- Cj_i) (68) 

Each of the n-1 points is represented this way, with j being the current point and j-1 the previous point. The 
first point j = 1 is the initial vehicle condition. The problem can be rearranged (Eq. (69)), so that it is a 
function of the previous point and the coefficient matrices. 


Xj = [I — 0.5 dtAj] 1 [(/ + O.bdtAj-i) Xj - i + 0.5 dt ( BjUj A Bj-iUj^.% + cj A Cj-i)\ (69) 

I is the identity matrix. A, B , and c are evaluated based on the previous iteration, prior to solving the 
optimization problem. This form replaces the continuous dynamics in Problem P3 as n-1 additional equality 
constraints levied on the problem. An advantage of formulating the equations in this manner is that the 
optimizer does not require a propagator. 

In the analysis a time step of 2.0 seconds was used. This was compared to smaller time steps of 1.0 
second and 0.5 second, with negligible change in the results. The smaller time steps increase the number of 
equality constraints levied on the problem, thus increasing the number of calculations and the run time of the 
simulation. Care should be taken to balance the accuracy from the propagation and the number of equality 
constraints. 


ANALYSIS AND RESULTS 


The entire problem was coded in Matlab to examine feasibility, convergence and descent trajectory trends. 
The convex optimizer chosen as the testbed was CVX. 19 CVX is a publicly available Matlab based solver for 
disciplined convex problems. The discrete version of Problem P3 is in a form accepted by CVX. Iterations 
were considered completed when the position vector between the previous and current iteration was less than 
3 cm apart at any trajectory point. 

The sample spacecraft used in the analysis consisted of 1000 kg dry mass with 400 kg of propellant. The 
trajectories used a small fraction of the usable propellant, so this can be thought of as having more dry mass 
and payload on the vehicle. The maximum thrust was 80 N and the minimum thrust was 20 N with a 225 
second Isp. 


Three sizes of a generic triaxial asteroid were investigated, Al: 1000 x 500 x 250 m, A2: 750 x 500 x 
250 m and A3: 500 x 500 x 250 m. Each asteroid rotates around the largest moment of inertia (+z axis) with 
a constant spin rate, which is reasonable as the majority of asteroids rotate around their largest moment of 
inertia. 17 The gravitational coefficients for the 2x2 triaxial spherical model can be analytically calculated and 
are given in Eq. (70). 17 


6 ^ 2,0 

6 ^ 2,2 





(70) 


The size of the asteroid determines these coefficients with a corresponding to the largest semi-major axis, 
/ 3 to the intermediate semi-major axis and 7 to the smallest semi-major axis. In addition, four different spin 
rates were examined. These are described by the amount of time it takes for the asteroid to complete one 
revolution. The periods are 12 hours, 8 hours, 4 hours and 2 hours. The analysis focused on the 8 hr period 
with the Al asteroid. The additional cases were confirmation of trends. 
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Table 1. Spacecraft initial conditions. 



ro x (m) 

ro y (m) 

ro z (m) 

v 0x (m/s) 

vo y (m/s) 

v 0z (m/s) 

NP 

-50 

-50 

1000+7 

2 

1 

0 

EQ 

1000+a 

-50 

-50 

0 

2 

1 



Figure 3. Vehicle trajectory for NP (top) and EQ (right) for a 380 second flight time. 


Four trajectories were analyzed during the investigation to ensure that the methodology is robust. The 
first two cases start 1000 m in altitude, uprange and out of plane from the landing site, along with an initial 
velocity needing to be slowed. One landing site is on the +z axis or the north pole of the asteroid (NP) with 
an equatorial trajectory landing on the +x axis (EQ). Figure 3 shows these two trajectories for a 380 second 
flight time. The initial spacecraft state is located in Table 1. The second set of cases hover 1000 m above the 
landing site and rotate with the landing site. There is both a north pole (NP_hov) and an equatorial (EQ_hov) 
hovering case. An additional challenge for the equatorial cases is the landing site rotation inducing a small 
velocity relative to the ACF frame. The north pole trajectories do not have this challenge. Its challenge 
is a stronger gravitational acceleration and a steep change in the gravitational acceleration approaching the 
surface. 
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Figure 4. Thrust profiles for EQ 400 s (left), 525 s (middle), 600 s, (right) flight time, 
showing the three different classes of thrust profiles. 
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The thrust magnitude follows the traditional bang-bang profile common in optimal control problems. Over- 
all there are three profiles that are prevalent, maximum-minimum-maximum, maximum-minimum and min- 
imum. Shorter flight times have the maximum-minimum-maximum profiles and long flight times have min- 
imum thrust for the entire flight. The maximum-minimum occurs in between the two. Figure 4 contains an 
example of each of these three classes of thrust profiles. 


NP EQ NP_hov EQ_hov 
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Figure 5. Propellant used as a function of flight time for 1000 x 500 x 250 m with an 8 hr period. 


Table 2. Flight times corresponding to events . 


(s) 

NP 

EQ 

NP_hov 

EQ_hov 

min prop 

488 

512 

501 

525 

start max-min 

465 

521 

474 

529 

start all min 

553 

598 

502 

536 


Table 3. Percent difference in mass from minimum propellant case. 



NP 

EQ 

NP_hov 

EQ_hov 

start max-min 

0.33% 

0.019% 

0.45% 

0.002% 

start all min 

1.8% 

2.5% 

0.17% 

0.19% 


Parameter sweeps changing the length of flight were performed for all the trajectories, asteroids and peri- 
ods. The parameter sweeps for the four trajectories landing on the 1000 x 500 x 250 m asteroid with an 8 hr 
period are shown in Figure 5. The propellant used is a convex function in the form of a quadratic spliced with 
a linear function. The point where the splice occurs (slope changes) and all periods use the same amount of 
propellant is where the thrust profile switches to minimum thrust for the entire trajectory. After this point the 
length of flight overwhelms the trajectory, gravity, asteroid size and period aspects in determining propellant 
usage. Table 2 lists the time that the switch from maximum-minimum-maximum to maximum-minimum 
occurs corresponding to this parameter sweep, along with the switch to all minimum thrust. 

Ideally, the spacecraft would fly the trajectory pertaining to the minimum propellant used over all flight 
times (optimal flight time). Table 2 lists the optimal flight time in the first row. This minimum propellant 
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trajectory occurs before the switch to all minimum thrust trajectories. For the EQ trajectory the minimum 
propellant case occurs before the switch to maximum-minimum, while it occurs after this switch for the 
other three trajectories. The propellent usage is flat near the minimum propellant trajectory, especially before 
the minimum propellant case. There is less than a 1% change in propellant used if the bum is 30 seconds 
shorter than the optimal flight time. After the burn, the propellant usage increases faster than before due to 
the proximity of the switch to the all minimum thmst on the hover cases. However, it remains less than a 
2% change within 10 seconds after the minimum propellant case. The switch to maximum-minimum occurs 
within 0.5% difference of the minimum propellant case as seen in Table 3. 



Figure 6. EQ trajectory 512 s flight time, for 1000 x 500 x 250 m with 8 hr period. 
Spacecraft position (top row) and velocity (middle row) relative to landing site. Thrust 
profile bottom row. 
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Figure 7. EQ 3-D position for 512 second flight time, for 1000 x 500 x 250 m with 8 hr period. 


A look at the fully designed trajectory for the minimum propellant EQ case for the 8 hr period 1000 x 500 
x 250 m can be found in Figure 6. The position and velocity are plotted with respect to landing site. The 
position path that the spacecraft will follow down to the surface is shown in Figure 7. This trajectory has zero 
velocity relative to the asteroid surface upon touchdown as was designed in the optimization problem. 
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Figure 8. Propellant used as a function of flight time for EQ 1000 x 500 x 250 m 
asteroid for all four periods. 
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Figure 9. Propellant used as a function of flight time for NP 1000 x 500 x 250 m 
asteroid for all four periods. 
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Figure 10. Propellant used as a function of flight time for EQ 2 hour flight period for 
the three asteroid sizes. 
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Figure 11. Propellant used as a function of flight time for NP 2 hour flight period for 
the three asteroid sizes. 


The trajectory design worked great for the three asteroid sizes and the four rotation periods for all four 
trajectories. The equatorial trajectories are sensitive to the rotation rate, with more propellant needed for the 
faster rotation (2 hours). Figure 8 shows EQ flight time parameter sweeps for the four revolution periods with 
the 1000 x 500 x 250 m asteroid. The other two asteroid sizes showed the same trends, as well as the EQ_hov 
trajectories. The north pole trajectories are not affected by the revolution speed of the asteroid, as the north 
pole landing site is on the rotation axis. Figure 9 confirms this, showing the parameter sweep for the 1000 x 
500 x 250 m asteroid for all four periods. Examining the trends for asteroid size, shows that the equatorial 
trajectories are not sensitive to the asteroid size, while the north pole trajectories are. Figure 10 depicts the 
EQ parameter sweep for a 2 hr period for the three asteroids, with very little difference between the three. 
The 8 hr period had no difference between the asteroids. Figure 1 1 contains the NP parameter sweep for a 2 
hr period for the three asteroids. As the asteroid decreases in size, the propellant usage increases. The hover 
trajectories showed identical trends to the non-hover trajectories. 

A good measure of how rapidly this methodology can work is how many iterations (successive solution 
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problems) were required to design the trajectory. Overall, between 2 to 5 iterations were required, with the 
majority being 3 or 4 iterations. Unlike propellant used, there were no trends among asteroid size, rotation 
speed or flight time that drove the number of iterations. The equatorial cases tended to have more trajectories 
that took 2 iterations and fewer 5 iterations. To give an idea of the simulation run time, a 600 second case with 
four iterations took 3.5 minutes, with each iteration solving for approximately 3300 variables. This code has 
not been scrubbed for efficiency, nor written in the most efficient language. It is expected that this time will 
decrease significantly, when programmed for efficiency. Also, in-depth analysis could reduce the number of 
nodes. 

CONCLUSION 

This paper demonstrates that the minimum fuel powered descent problem can be turned into a convex 
optimization problem, while still retaining the same optimal solution as the original problem. Successive 
solution method is the key to handling the nonlinear gravity model and adds the flexibility of switching to 
higher fidelity gravity models in the future. The algorithm successfully handled a wide range of trajectories 
fleshing out its ability to be adaptable to the asteroid parameters after the asteroid has been characterized upon 
spacecraft arrival. Trajectories landing near the asteroid’s equator are sensitive to the rotation speed, while 
trajectories landing near the poles will be sensitive to the asteroid size in terms of propellant required for 
landing. Propellant used for flight times near the minimum propellant case differ very little from the optimal 
case. Therefore, the exact flight time of the minimum propellant case does not need to be known with a 
high degree of accuracy. Overall, this is a viable algorithm for rapidly designing powered descent trajectories 
autonomously on-board the spacecraft. 
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